Random transition-rate matrices for the master equation 
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Random-matrix theory is applied to transition-rate matrices in the Pauli master equation. We 
study the distribution and correlations of eigenvalues, which govern the dynamics of complex stochas- 
tic systems. Both the cases of identical and of independent rates of forward and backward transitions 
are considered. The first case leads to symmetric transition-rate matrices, whereas the second cor- 
responds to general, asymmetric matrices. The resulting matrix ensembles are different from the 
standard ensembles and show different eigenvalue distributions. For example, the fraction of real 
eigenvalues scales anomalously with matrix dimension in the asymmetric case. 
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I. INTRODUCTION 

The Pauli master equation is encountered in many 
fields of science such as physics, chemistry, and biology. 
It describes the time evolution of probabilites for a sys- 
tem to be in certain states. Formally identical rate equa- 
tions describe the dynamics of concentrations or popula- 
tions of certain entities. The dynamics of probabilities is 
described by the Pauli master equation 



Pi 



(1) 



where P, is the probability to find the system in state 
i = 1, . . . , N and Rij is the transition rate from state j 
to state i. Evidently, the rates of change of probabilities 
depend only on the probabilities at time t, i.e., Eq. (Q]) 
describes a memory-less or Markovian process. Equation 
(P) ensures that the total probability is conserved, 



RjiPi) = 0. 



(2) 



Typical applications in physics include lasers hj, disor- 
dered conductors microelectronic devices [3[, quan- 
tum dots [4[ , and molecular electronics [5] . In these cases 
one can, in principle, obtain the Pauli master equation by 
first deriving a quantum master equation for the reduced 
density matrix of a small system, which is obtained by 
tracing out the reservoir de gree s of freedom from the full 
density operator [|| 0, H, [1,113] ■ If the off-diagonal com- 
ponents of the reduced density matrix decay rapidly, it 
is sufficient to keep only the diagonal components rep- 
resenting the probabilities Pi of states \i) of the small 
system. In certain fields, for example in transport and 
laser theory, the resulting Eqs. ([I]) are often called rate 
equations. 

However, if even the small system is complicated, such 
as a system of interacting enzymes, this route becomes 



unfeasible. In applications outside of physics, where i 
could refer to the state of a technical or social process, 
a quantum-statistical description becomes inappropriate 
in any case. One would then view Eq. (TT]) as the funda- 
mental description. 

Our goal is to make progress in the understanding of 
the master equations for complex systems. The number 
N of possible states will typically be large. It should be 
noted however that complex behavior can already emerge 
for moderate N. An example is provided by the differ- 
ential conductance calculated in Ref. [llj for a magnetic 
molecule with magnetic anisotropy axis not aligned with 
the applied magnetic field, where N = 20, but due to 
noncommuting terms in the Hamiltonian many rates are 
nonzero and are distributed over a broad range. 



A. Properties of the master equation 

We first recount some basic properties. It is clear that 
one can rewrite Eq. |T]) in the form 



(3) 



or P = AP with the transition-rate matrix, or, for short, 
rate matrix, 



A, 



R, 



for i =/= j 



It follows that the column sums vanish, 
y^Aij = for all j. 



(4) 



(5) 



Note that (d/dt) J2i Pi — J2ij AijPj vanishes for all Pj 
if and only if Eq. |(SJ| holds. The constraint (0 is thus 
dictated by conservation of probability. From Eq. Q it 
is also clear that 
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Aij > for all i ^ j 



(6) 
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if we interpret the Rij as transition rates. A matrix sat- 
isfying the inequalities ^ and Y^i Aij — f° r au j i s 
called a compartmental matrix. 



Equation ([3]) can be solved by the ansatz P 



e At v, 

which leads to the eigenvalue equation Av — Av. Since 
A is generally not symmetric, the eigenvalues A and the 
components of the right eigenvectors v can be complex. 
However, since A is real, the equation Av = Av implies 
Av* = A*v*. Thus, the eigenvalues are real with real 
eigenvectors or form complex conjugate pairs with their 
eigenvectors also being complex conjugates. 

Let v„ be the right eigenvector to eigenvalue A„. It 
is well known that there is always at least one strictly 
zero eigenvalue, which we call Ao = 0: the constraint (JSJ) 
implies that A has a left eigenvector (1, 1, . . . , 1) to the 
eigenvalue Ao = 0. The corresponding right eigenvector 
Vo describes the stationary state. 

A real eigenvector v„ with real eigenvalue A„ describes 
a contribution to the probability vector P that decays 
exponentially with the rate — A n . A complex conjugate 
pair of eigenvectors v„, v* with eigenvalues A„, A* can 
be combined to form the two independent real solutions 



(e A "*v„ + e A "*v;)/2 and (e A "*v„ - e A: * 
the components of v„ as v n j — v^e 1 ^" 
obtain the solutions 



v*)/2i. Writing 
with v„j real, we 



ReA„t 

v n] e x 



cos(Im X n t 
sin(Im X n t - 



The initial values at time t = are clearly Re v n j and 
lmv n j, respectively. We thus find damped harmonic 
oscillations with damping rate —Re A„ and angular fre- 
quency ImA„. We obtain the solution at all times by ex- 
panding the initial probability vector P(t = 0) into the 
basis of real vectors v„ (for real A n ) and Rev„, Imv„ 
(for complex conjugate pairs A„, A*). 

An eigenvalue A„ with Re A„ > would be unphysical, 
since the corresponding contribution to the probabilities 
would diverge for t — > oo. However, for any compartmen- 
tal matrix the spectrum is contained in {A|Re A < 0}u{0} 
[T2I [T3I ] . Thus all eigenvalues are either zero or have a 
strictly negative real part. 

The Perron-Frobenius theorem [3, EH] applied to the 
non-negative matrix A — a mm I, where a m - m < is the 
minimum of An and I is the N x N unit matrix, shows 
that the right eigenvector Vo to Ao has only non-negative 
components. This ensures that the probabilities in the 
stationary state are non-negative. 



B. Random rate matrices 

As noted above, even relatively simple problems lead to 
master equations with rates Aij, i ^ j, distributed over a 
broad range. In problems with large numbers of states it 
is often impractical to obtain all independent components 
A^. This situation is reminiscent of Hamiltonians for 
complex systems. Difficult problems of this type concern 



atomic nuclei and quantum dots, where the Hamiltonian 
is too complicated to write down explicitly, but cannot 
be simplified by methods restricted to weakly interact- 
ing systems. For these systems, random-matrix theory 
(RMT) 0, El El, EH has lead to significant progress. 
The main assumption is that a Hamiltonian of this type 
is a typical representative of an ensemble of Hamiltoni- 
ans of appropriate symmetry. While this approach does 
not allow one to obtain specific eigenvalues, it does pro- 
vide information about the statistical properties of the 
spectrum 0, El El E! • 

Our point of departure is to treat the rate matrix A 
for a complex system as an element of a suitable random- 
matrix ensemble. In the case of transport through quan- 
tum dots, this is complementary to treating the Hamil- 
tonian of the quantum dot as a random matrix, which 
has been done extensively [l7j . 

Since the rate matrix A must satisfy the conditions 
([5]) and © , we define the exponential general rate-matrix 
ensemble (EGRE): The EGRE is formed by real N x N 
matrices A with independently identically distributed off- 
diagonal components A^ with the distribution function 



p(A ij ) = { (R) e h»A,>« 
otherwise 



CO and the diagonal components 



(8) 



(9) 



The exponential distribution of rates A^ is viewed as the 
least biased distribution of non-negative numbers. We 
will also present results that do not depend on the spe- 
cific distribution function p. We will see that the spe- 
cific distribution becomes irrelevant in the limit of large 
N, at least if all its moments exist. The distribution of 
components is thus not the most fundamental difference 
between the EGRE and the well-known ensembles stud- 
ied in the context of random Hamiltonians. Rather, one 
such difference lies in the constraint {5| or ^ . The other 
is that the rate matrices are real but not symmetric and 
thus not hermitian [20| . 

Ensembles of non-hermitian matrices have been stud- 
ied in detail, starting with Ginibre's work on Gaussian 
ensembles of non-hermitian matrices with real, complex, 
and quaternion components [21| . We will compare our 
results to the real Ginibre ensemble. 

To be able to analyze the importance of the asymme- 
try, we also define the exponential symmetric rate-matrix 
ensemble (ESRE) : The ESRE is formed by real symmet- 
ric N x N matrices A with independently identically dis- 
tributed components A^ above the diagonal (i < j) with 
the distribution function given by Eq. ([5]) and the diag- 
onal components given by Eq. ([9]). 

Another possible choice is a two- valued distribution of 
rates, where a transition from state j to state i is either 
possible or impossible, and all possible transitions have 
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the same rate. This case with symmetric rates has been 
studied by various authors [H, l23l [24| . It is essentially 
equivalent to adjacency matrices of random simple net- 
works. 

An ensemble of real symmetric matrices satisfying 
Eq. ([5]) but with a Gaussian distribution of Aij has also 
been studied [24[ . This case cannot easily be interpreted 
in terms of a master equation, since the can be neg- 
ative. We will compare our results for the eigenvalue 
spectrum to these works below. 

The remainder of this paper is organized as follows: In 
Sec. |TT] we consider the simpler case of symmetric rate 
matrices (the ESRE) and obtain results for the eigen- 
value density and for the correlations between neighbor- 
ing eigenvalues. In Sec. IIIII we then study general rate 
matrices (the EGRE) and obtain results for the eigen- 
value density, now in the complex plane, and for the 
correlations of neighboring eigenvalues. We conclude in 
Sec. lIVl A number of analytical derivations are relegated 
to appendices. 



II. SYMMETRIC RATE-MATRIX ENSEMBLE 

We first consider ensembles of symmetric rate matrices 
A. These describe processes where transitions from any 
state j to state i and from i to j occur with the same 
rate, A^j = Aj^. 



A. Spectrum 

As noted above, the spectrum always contains the 
eigenvalue Ao = 0. The corresponding eigenvector for 
symmetric matrices is (1,1,..., 1) or, normalized to unit 
probability, (1/N,1/N, . . . ,1/N). For symmetric rates, 
the stationary state is thus characterized by equal distri- 
bution over all states i. We are interested in the distribu- 
tion of the other eigenvalues A n , n = 1, . . . , N — 1, which 
are all real. We have also seen in Sec. II Al that A„ < 0. 
Since there is no further constraint, the probability of A„ 
for any n > being exactly zero vanishes. 

To simplify the calculations, we shift the matrices so 
that they have zero mean. We discuss this immediately 
for general matrices. Also, nothing here depends on the 
distribution function p of the rates Aij, as long as the 
average (R) = (A^) exists. We define 



the eigenvalue Ao = 0, we have wj^v n = 0. Since 



A = A-(A), 



(10) 



where here and in the following angular brackets de- 
note the average over the matrix ensemble under con- 
sideration. Here, (A) has the components (Aij) = (R) 
for i ^ j and (An) = -(N - 1) (R). Is follows that 
J2i Aij = for all j. Consequently, A has a left eigen- 
vector -Wq = (1, 1, . . . , 1) to the eigenvalue Ao = 0. 

Let v n be the right eigenvectors of A to the eigenvalues 
X n , n = 1, . . . , N — 1. Since is the left eigenvector to 




v„ is a right eigenvector of (A) to the eigenvalue —N(R). 
Therefore, v„ is also a right eigenvector of A to the eigen- 
value A„ = A„ + N(R). The result is that the shifted 
matrices A also have one eigenvalue Ao = and that 
the remaining eigenvalues are just the eigenvalues of A, 
shifted by N(R). 

We now derive the average of eigenvalues A n , here and 
in the following excluding Ao — 0. We have (A)' = (A)' — 
N(R), where angular brackets with a prime denote the 
average over all eigenvalues, excluding the exact zero. 
Since this leaves N — 1 eigenvalues, their average is the 
trace of the matrix, to which the zero eigenvalue does not 
contribute, divided by N — 1. Consequently, 



1 



1 



<A)' = ] ^ T Tr(A> = F - T Tr0 = (12) 



so that 



(A)' = -N(R). 



(13) 



This result is independent of the specific distribution 
function of rates, p, as long as (R) exists. 

We next calculate the low-order central moments 

Mm = (\ m Y = ((A - (A)')" 1 )' = ((A + N(R)) m Y (14) 

of the eigenvalues A„, n > 0. The central moments are 
identical to the central moments of the shifted values X n . 
Unless otherwise noted, our results for fi m hold for an 
arbitrary distribution function of rates, p, as long as the 
moments exist. It is instructive to show the calculation 
of the second moment explicitly. We find 



= (A 2 ) 



1 

N- 1 



Tr(l 2 



1 



N _ 1 Z_^\^*]~Ti 



N 



(AijAji) 



3+i 



Using Aij — Aji and J2k ^ki — 0, we obtain 



j^i E (e<4> + E 



(16) 



With (A^) = we finally get 



i jjti 



where (5R 2 ) = (A'fj) — (Aij) 2 for i ^ j is the second 
central moment of p(Aij). For the special case of an 
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exponential distribution we have (SR 2 ) = (R) 2 and thus 
p 2 = 2N(R) 2 . 

The important consequence is that while the mean of 
the nonzero eigenvalues of the unshifted matrices A scales 
with N, Eq. (|13p . the width of their distribution is only 
^ = y/2N (SR 2 ) oc VN. Thus for large N the distri- 
bution of eigenvalues contains a single eigenvalue Ao = 
and the remaining N — 1 eigenvalues form a narrow dis- 
tribution around —N(R). In physical terms, nearly all 
deviations from the stationary state decay on the same 
time scale 1/N(R). 

All moments can be obtained by the same method: We 
first write the average in terms of a trace, split the sum 
into terms with equal or distinct matrix indices, and use 
^2 k Aki = 0. With Aij = Aji and (Ay) = we obtain 
the moments. Since the enumeration of all possible cases 
of equal or distinct indices is cumbersome, we have used a 
symbolic algebra scheme implemented with Mathematica 
[25| . The results up to to = 8 are shown in Table [T] for 
a general distribution. The moments are expressed in 
terms of the central moments (5R n ) = ((Aij — (Aij)) n ). 
Note that in the limit of large N, the moments fi m for 
even to only depend on the second moment (5R 2 ). We 
will return to this point shortly. 

Table [TT1 shows the central moments fx m up to to = 10 
for the exponential distribution of Ay, i < j (ESRE). For 
the exponential distribution, one has (SR n ) — In (R) n , 
where In = n\ Efc=o(~l) fc /^! * s ^ ne subfactorial. Table 
HJalso contains the leading large- N terms for the ESRE. 
At least up to to — 10, the even moments scale as u m ~ 
£ or i ar g e jy^ as expected from the scaling of p, 2 - 
However, the odd moments scale only as p, m ~ 7V( m—1 " 2 . 
If this holds for all to, the distribution of A approaches 
an even function for large TV. This is indeed the case, as 
we shall see. 

The density of eigenvalues A„ can be obtained from the 
resolvent [26| G(z) = (z — A)^ 1 . The density is given by 
the spectral function 



Pall (z) 



1 

'tvN 



ImTr (G(z + in)}, 



(18) 



where n — ► + at the end of the calculation. The density 
includes the exact zero eigenvalue so that we can write 



Pall (z) 



1 

N 



S(z) 



N -1 

N 



P(z), 



(19) 



where p(z) is the normalized density of nonzero eigen- 
values. In the limit of large N, the eigenvalue density 
p & w(z) = p(z) only depends on the second moment (SR 2 ) 
of the distribution function p of rates, at least as long as 
all moments of p exist. The proof is sketched in App. [A] 
That the eigenvalue distribution generically becomes in- 
dependent of p for large N has been conjectured by Mehta 
(conjecture 1.2.1 in Ref. [lj|). However, the second part 
of this conjecture, stating that the density of eigenval- 
ues is the same as for the Gaussian orthogonal ensemble 
(GOE), is not true for our ensemble. 



Since the density of eigenvalues A„, n > 0, of the 
shifted matrices A only depends on the second moment 
(SR 2 ) for large N, we can obtain the large-TV behav- 
ior from any distribution with that second moment. We 
choose the Gaussian distribution 



Pg(A, 



1 



VMSR 2 ) 



exp 




(20) 



For this distribution together with the constraint 
Y^j AjA — 0, the eigenvalue density is known for large 
TV |24f : The averaged resolvent is the solution of 



(G(z)) 



where 



y/N(SR 



9(z) 



1 g z- N (SR 2 )(G( z) )\ 



i 



dx ■ 



-x 2 /2 



VZirJ-co z-x 
This integral can be evaluated, 



(22) 



g{z) 



i 



-2 + erfc • 



V2 ; 

g(z) has a cut along the whole real axis. The density 
p(z) is thus nonzero for all real z. Equations (fT5| and 
(l2"Tj) imply that yj N(5R 2 ) p(z) is a universal function of 
z j yjN (SR 2 ). The same distribution in the large- N limit 
was found for adjacency matrices [22j, |23|. The corre- 
spondingresult for the GOE is the well-known semicircle 
law [l6l Tl9l]. It is worth pointing out that the differ- 
ent eigenvalue density results only from the constraint 

E» hi = °- 

We now study the eigenvalue density for the ESRE 
for finite N. We perform Monte Carlo simulations by 
generating a number n r of realizations of matrices from 
the ESRE for given N, shifted according to Eq. (JTUJ). 
The matrices are diagonalized and the eigenvalue with 
the numerically smallest magnitude, which corresponds 
to Ao = 0, is droppe d. The e igenvalues are rescaled ac- 
cording to A — > \/y / N(SR 2 ). Finally, histo grams with 
500 bins are generated. 

Results for N = 2, 10, 100, 1000, 10000, and 00 are 
shown in Fig. [TJ For N — > 00, we solve Eq. (|2"Tj) . For 
N = 2, the matrices have a single nonzero eigenvalue 
— 2A12 with distribution following from Eq. ©. For each 
of the other values of N, n r N — 10 7 eigenvalues have 
been generated. Figure [T] shows that the distribution 
changes smoothly from shifted exponential for N = 2 
to the known universal function for N — > 00. The inset 
in Fig. [1] shows the unsealed eigenvalue density of the 
unshifted ESRE to illustrate that the mean scales with 
N, whereas the width scales with s/N. 

While we have shown that nearly all nonzero eigenval- 
ues lie in a narrow interval around their mean for large 
TV, the dynamics after a transient will be dominated by 
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TABLE I: Central moments fi m , m = 2, ...,8, of the nonzero eigenvalues A for ensembles of symmetric rate matrices. The 
results hold independently of the distribution function p of rates Aij, i < j, as long as the moments exist. Here, (5R n ) is the 
n-th central moment of p. 



m fi m (symmetric matrices, general distribution) 

2 2N{SR 2 ) 

3 -4TV(5i? 3 ) 

4 7V[9(7V -2)(SR 2 ) 2 + 8(SR 4 }] 

5 -2N[25(N -2)(SR 2 ){SR 3 ) + 8(5R 5 )] 

6 TV[4(14TV 2 - 737V + 90)(<57? 2 ) 3 + 73(7V - 2)(SR 3 } 2 + 132(7V - 2)(5R 2 )(SR 4 } + 32(6R 6 }] 

7 -27V[7(417V 2 - 2117V + 258) (5 R 2 ) 2 (SR 3 ) + 203(7V - 2) (S R 3 } {5 R 4 } + 168(7V - 2)(SR 2 )(SR 5 } + 32(SR 7 }] 

8 7V[(4317V 3 - 40427V 2 + 120217V - 11322){<5i? 2 } 4 + 6(3067V 2 - 15617V + 1898)(5 J R 2 ) 2 { ( S J R 4 ) + 593(7V - 2){8R 4 ) 2 

+ 1088(7V - 2)(SR 3 )(5R 5 ) + 4(7V - 2)(5077V - 1574)(<5i? 2 )(<5i? 3 ) 2 + 832(7V - 2)(5R 2 ){5R 6 ) + 128(5R 8 )} 



TABLE II: Second column: central moments p m , m — 2, ... ,10, of the nonzero eigenvalues A for ensembles of symmetric rate 
matrices, assuming an exponential distribution of rates (ESRE). Third column: leading term of fi m for large TV. 



m 


Mm (ESRE) 


fi m (ESRE, TV > 1) 


2 


27V(i?) 2 


2N{R) 2 


3 


-87V(i?) 3 


~8N{R) 3 


4 


97V(7V + 6)(.R) 4 


97V 2 (i?> 4 


5 


-47V(257V + 126)(R) 5 


-1007V 2 (i?) 5 


6 


47V(147V 2 + 297TV + 1470) (J?) 6 


567V 3 (R) 6 


7 


-47V(2877V 2 + 40467V + 20424) (R) 7 


-11487V 3 (R) 7 


8 


7V(4317V 3 + 205947V 2 + 2505767V + 1311648) (R) 8 


4317V 4 (Rf 


9 


-47V(34537V 3 + 950217V 2 + 10894147V + 5957208) (J?) 9 


-138127V 4 (R) 9 


10 


27V(19717V 4 + 1726 577V 3 + 37371277V 2 + 421066107V + 241175496) (J?) 10 


39427V 5 (R) 10 



the slowest process. The slowest non-stationary process 
is governed by the eigenvalue Ai < which is smallest 
in magnitude. It is conceivable that matrices from the 
ESRE typically have an eigenvalue Ai close to zero. For 
example, Ai could scale with a lower power of TV com- 
pared to the mean —N(R) . If the fraction of such anoma- 
lously slow rates decreased for large TV, they might not 
be visible in the density plots in Fig. [T] 

To check this, we plot the mean (Ai) as a function of 
TV in Fig. [2] The average slowest rate |(Ai)| is signifi- 
cantly smaller than the average rate |(A)'| for small TV, 
as one would expect from the width y/Ji2 oc y/~N. On the 
other hand, for large TV, |(Ai)| approaches |(A)'|. Thus 
we do not find evidence for anomalously slow processes. 
Instead, the slowest rate is consistent with the mean and 
width of the eigenvalue distribution p(A). 

B. Eigenvalue correlations 

Since the eigenvalue density for the ESRE differs sig- 
nificantly from the GOE, one might ask whether the cor- 
relations between eigenvalues are also different. In the 
GOE, the distribution function of differences of neigh- 
boring eigenvalues A, A' approaches zero as |A' — A| for 



Figure [3] shows the distribution function pnn(AA) of 
separations AA = A n +i — A„ of neighboring eigenval- 
ues for the ESRE (here, the A n are assumed to be or- 
dered by value) . The zero eigenvalue Ao = is excluded. 
Since the width of the eigenvalue distribution scales as 
V~N, while the number of eigenvalues for a given realiza- 
tion scales as TV, the typical separati on should scale as 
1/y/N. We therefore rescale AA -> y / N/(R} 2 AA. Fig- 
ure [3] shows that the rescaled distribution approaches a 
limiting form for TV — * oo. Furthermore, the distribution 
function pnn(AA) is linear in AA for small AA for all TV. 
Thus the distribution of nearest-nei ghb or separations be- 
haves essentially like for the GOE The constraint 
(|5|), which is responsible for the deviation of the eigen- 
value distribution from the GOE result, does not have a 
comparably strong effect on the eigenvalue correlations. 
The reason is very likely that the joint probability distri- 
bution p(Xi, A2, . . • , Aat_i) of the eigenvalues [19(, while 
being complicated for the ESRE, does contain the factor 
n„n',o<n<n' l^n - A„' | , which determines the exponent 
(3 = 1 in pnn ~ AA* 3 . 
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FIG. 1: (Color online) Scaled density of nonzero eigenvalues 
of shifted symmetric rate matrices A = A — (A) . The results 
for N = 2 and N — * oo are exact, see text. The curves for 
N = 10, 100, 1000, 10000 are histograms with 500 bins for 
10 7 eigenvalues for matrices randomly chosen from the ESRE. 
Inset: unsealed distribution of eigenvalues of the unshifted 
matrices A for N = 2, 10, 100. 
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FIG. 3: (Color online) Scaled distribution of nearest-neighbor 
separations AA of nonzero eigenvalues for the ESRE for N = 
10, 100, 1000, 10000, from the same data sets as in Fig. Q] 
The curve for N = 1000 is nearly obscured by the one for 
N = 10000. 



A. Spectrum 
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FIG. 2: (Color online) Average smallest in magnitude eigen- 
value, (Ai), of matrices from the ESRE, as a function of 
N . The open circles denote numerical results for n r = 5000 
(1000) realizations for N < 1000 (N > 2000). Error bars de- 
noting the statistical errors are shown. The filled square de- 
notes the result Ai = — 2(R) for N — 2. The dashed straight 
line denotes the mean of nonzero eigenvalues, —N{R). 



III. GENERAL RATE-MATRIX ENSEMBLE 



We now turn to the ensemble of general, asymmetric 
rate matrices (EGRE). Compared to the ESRE, it de- 
scribes the opposite extreme of independent rates Aij 
and Aji for forward and backward transitions. 



As noted, there always exists an eigenvalue Ao = with 
left eigenvector (1, 1, . . . , 1) . Other than for the symmet- 
ric case, the corresponding right eigenvector is different. 
We are interested in the distribution of the other eigen- 
values X n , n = 1 , . . . , N — 1 , which are now complex with 
negative real parts. We have already shown in Sec. |TT] 
that the mean of nonzero eigenvalues equals —N(R), see 
Eq. (p~5|) . We shift the matrices according to Eq. (fTT)|) so 
that they have zero mean. 

We define the expectation values 



(A* 



((A -(A))' 



((A + N{R)) 7 



(24) 



in analogy to the ESRE, but they are not the central 
moments of the distribution of nonzero eigenvalues. In- 
stead, the central moments have to be defined for a two- 
dimensional distribution in the complex plane, 



^ mn = ((ReX + N(R)) m {ImXy 



(25) 



Since the eigenvalues are real or form complex conjugate 
pairs, we have [i mn — for odd n. We show in App. [B] 
that the shifted eigenvalue distribution only depends on 
the second moment (5R 2 ) of p, like we found for the 
symmetric case. We here call the ^ m in Eq. ([24]) the 
pseudomoments. They are all real, since the eigenvalues 
are real or form complex conjugate pairs. 

The pseudomoments p m can be obtained in the same 
way as for symmetric matrices. The results are different, 
since (AijAji) — (SR 2 ) for the symmetric case, whereas 
(AijAji) = for the general case. We present the pseu- 
domoments fi m up to m — 8 for a general distribution 
function p(Aij) in Table UlTl and up to m = 10 for the ex- 
ponential distribution (EGRE) in Table |Wj The scaling 
of fi m for even and odd m and large N is the same as for 
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the ESRE. In the limit N — > oo, only the even pseudo- 
moments survive. Interestingly, at least up to m = 10, 
these agree with the central moments of a real Gaus- 
sian distribution, /x^ = (to — 1)!! (N (8R 2 )) m / 2 , where 
nil = n(n — 2)(n — 4) . . . is the double factorial. We show 
in App. Othat this identity holds for all even to. 

The eigenvalue distribution in the complex plane can 
be obtained from the non-analyticities of the averaged 
resolvent (G(z)) = {{z - A)" 1 ) [H, [27j]. However, un- 
like for symmetric matrices, the non-analyticities are not 
limited to a branch cut along the real axis. For what 
follows, it is more convenient to employ the method of 
hcrmitization [13. We define the 2N x 2N matrix 



N = 20 



0.5 - 



A 
V 

^ 0.0 








-0.5 - 




H(z,z*) 





A T - z*I 



A- zl 



(26) 



-4 



-2 

Re X I (N<R> 2 



2 
1/2 



where A T is the transpose of A. H(z, z*) is hermitian for 
any complex z. With the resolvent of 7i, 



g(r};z,z*) = 



1 



rj — H(z, z*) ' 

the density of eigenvalues in the complex plane is \2 



(27) 



Pa n(x,y) = ^^Tr 2N r Q jWo;*,**)), (28) 



where z = x + iy, the derivative with respect to z* is to 
be taken with z fixed, and Tt 2 n denotes the trace over 
a 2N x 2N matrix. Using this representation, we show 
that for large N the eigenvalue density only depends on 
the second central moment (SR 2 ) of the distribution of 
rates Aij. The proof is sketched in App. [BJ Edelman et 
al. [29( have conjectured that this is generically the case 
for asymmetric matrices. 

We now present numerical results for p(x, y) for the 
EGRE, as a function of the matrix dimensions N. As 
above, p a n contains all eigenvalues, whereas p excludes 
the exact zero. We will compare the results to the Gini- 
bre ensemble of real asymmetric matrices with Gaussian 
distribution of components (Ginibre orthogonal ensem- 
ble, GinOE) [Hi HI, [H, Hfl, U H 13, which is the 
closest relative of the EGRE that has been studied in 
detail. 

As observed above, the eigenvalues A of A can be ei- 
ther real or form complex conjugate pairs. The numerical 
simulations show that both types of eigenvalues indeed 
occur. A typical eigenvalue density is shown in Fig. 0] 
for N — 20. We assume that the square root of the sec- 
ond pseudomoment, ^fpa — y / N(R) 2 , describes the typ- 
ical width of the distribution and rescale the eigenvalue 
density accordingly. The real and complex eigenvalues 
are clearly visible. Here and in the following "complex" 
should be understood as "not real." Figure [4] already 
suggests that the distribution of nonzero eigenvalues of 
A becomes a narrow peak around —N(R) for large N, 
like for the ESRE. We return to this point below. 



FIG. 4: (Color online) Scaled distribution function of nonzero 
eigenvalues A of shifted general rate matrices A of dimension 
N = 20. More specifically, a two-dimensional histogram with 
500 x 500 bins was populated for n r matrices randomly chosen 
from the EGRE, with n r N = 4 x 10 7 . 



The question arises of what fraction /r of the nonzero 
eigenvalues are real. For the GinOE, this fraction is 
known analytically 29]. (The probability of finding ex- 
actly Nth real eigenvalues for N x TV matrices from the 
GinOE is also known [32].) Edelman et al. [2§] de- 
rive various equivalent expressions for the expected num- 
ber of real eigenvalues, (JVr), from which we obtain 
jGmOE _ ^tVr)/^. We here quote an expression in terms 
of the hypergeometric function 2 F\ 29]: 



GinOE 



1 

2N 



2 T{N + 1/2) 



7T r(jv + i) 

For large N, this becomes (29[ 



•Fx 1, 



2' ' 2 



fi 



GinOE 



2 



(29) 



(30) 



For the GinOE, the fraction of real eigenvalues thus 
asymptotically decays with a simple exponent of —1/2. 

Figure [5] shows the fraction /r as a function of N for 
the EGRE and for comparison the exact result for the 
GinOE. For N = 2, /r must be unity, since the single 
nonzero eigenvalue cannot be a complex conjugate pair. 
The results clearly differ from the GinOE and decay more 
slowly for large N. A fit of a power law /r ~ foN~ a to 
the data points for N = 2000 and 5000 is also included in 
Fig. El We obtain / » 1.37 and a w 0.460. The large-N 
behavior is inconsistent with the exponent 1/2 found for 
the GinOE. This is remarkable, since all other scaling re- 
lations we have so far found, as well as the ones for the 
GinOE, only contain integer powers of >/~N. Physically, 
this means that the fraction of eigenvectors describing 
purely exponentially decaying deviations from the sta- 
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TABLE III: Pseudomoments (x m , m = 2, . . . , 8, of the nonzero eigenvalues A for ensembles of general rate matrices. The results 
hold independently of the distribution function p of rates Aij, i 7^ j, as long as the moments exist. 



m fj, m (general matrices, general distribution) 

~2 N{5R 2 } 

3 —N(5R 3 ) 

4 N[3{N -1)(8R 2 ) 2 + (5R 4 )] 

5 -7V[10(7V - 1)(SR 2 )(SR 3 ) + (SR 5 )] 

6 7V[(157V 2 - 497V + 38)(5R 2 ) 3 + 10(7V - 1)(SR 3 ) 2 + 15(7V - 1)(SR 2 )(SR 4 ) + (SR 6 )] 

7 -7V[21(57V 2 - 177V + 14) (SR 2 ) 2 (5 R 3 ) + 35(iV - 1)(SR 3 )(5R 4 ) + 21 (TV - 1)(SR 2 )(8R 5 ) + (SR 7 )] 

8 7V[3(357V 3 - 2407V 2 + 5517V - 422)(SR 2 ) 4 + 6(357V 2 - 1217V + 102) (SR 2 ) 2 (SR 4 ) + 35(7V - 1)(5R 4 ) 2 

+56(7V - 1)(SR S ){SR 5 ) + 56(57V 2 - 187V + 16) (8 R 2 )(SR 3 ) 2 + 28(7V - 1)(SR 2 )(SR 6 ) + (SR S )] 



TABLE IV: Second column: pseudomoments fi m , m = 2, . . . , 10, of the nonzero eigenvalues A for ensembles of general rate 
matrices, assuming an exponential distribution of rates (EGRE). Third column: leading term of fi m for large TV. 



m 


lim (EGRE) 


ix m (EGRE, TV > 1) 


2 


N(R) 2 


7V(7?) 2 


3 


~2N(R) 3 


-27V(i?) 3 


4 


3N(N + 2)(R) 4 


37V 2 (R) 4 


5 


-4N(5N + 6)<i?) 5 


-207V 2 (R) 5 


6 


7V(15iV 2 + 126/V + 128)(ii> 6 


157V 3 (R) 6 


7 


-67V(35/V" 2 + 1407V + 148) (R) 7 


-2WN 3 {R) 7 


8 


7V(1057V 3 + 22907V 2 + 62707V + 7476)(i?) 8 


1057V 4 (R) 8 


9 


-87V(3157V 3 + 29537V 2 + 67417V + 9018) (R) 9 


-25207V 4 (R) 9 


10 


7V(9457V 4 + 424947V 3 + 2491747V 2 + 5328407V + 774744) (R) 10 


9457V 5 (R) 10 



tionary state scales with a nontrivial power —a of the 
number of states. 

To pinpoint the origin of the anomalous scaling, we 
have also evaluated /r for ensembles of matrices of di- 
mension TV = 5, 50, 500 satisfying the constraint ([5|), but 
with Gaussian distribution of rates Aij, i 7^ j. This is the 
asymmetric analogue of the symmetric ensemble studied 
by Staring et al. [24]]. The results are shown as crosses 
in Fig. [5] They clearly approach the EGRE results for 
large TV, not the GinOE. It is thus the constraint ([5]) that 
leads to the anomalous scaling. 

In the following, we will consider the real and com- 
plex eigenvalues separately. Figure [6] shows the den- 
sity /Oh of shifted real nonzero eigenvalues A, normal- 
ized to unity and rescaled with the square root of the 
pseudomoment ^/JH = V / N(R) 2 , for TV = 2, 10, 100, 
1000, 5000. For TV = 2, the single nonzero eigenvalue is 
A = — A12 — A2i- In the EGRE, its distribution function 
is Pk(A) = (2/(i?) - \/(R) 2 ) exp(A/(i?) -2) for A < 2(R) 
and zero otherwise. For the other values of TV, Fig. [S] 
shows numerical results. The noise increases for large TV, 
not only because n r 7V was smaller for TV = 5000 but also 
because /r decreases with increasing TV. It is obvious 
however that the distribution for large TV is quite differ- 
ent from the eigenvalue density for the ESRE, Fig. Q] 

The distribution clearly becomes more symmetric for 
TV — > 00, as it must, since the large- TV result only de- 



pends on the width of the distribution of rates Aij. 
There is an indication that the distribution develops non- 
analyticities with sudden changes of slope in the limit 
TV — > 00. This is not unexpected, since the scaled distri- 
bution of real eigenvalues of the GinOE is uniform on the 
interval [—1,1] and zero otherwise [29|, [33| and thus also 
shows non-analyticities. Compared to the ESRE (Fig.[T]), 
the convergence to the large-TV limit is slower for the 
EGRE (Fig. [5]). In fact, from Fig. [S] we cannot exclude 
the possibility that the width scales with an anomalous 
power of TV, different from 1/2. 

Turning to complex eigenvalues, we note that for large 
TV nearly all eigenvalues belong to this class, since the 
fraction /r of real eigenvalues approaches zero. We plot 
their distribution function pc in the complex plane for 
TV = 100 and 2000 in Fig. The scaled distribution for 
TV = 5000 is virtually indistinguishable from the one for 
TV = 2000. From Figs. Hand [3 we see that the distribu- 
tion becomes more symmetric with respect to inversion 
of the real part as TV increases. 

The widths of the distribution in the real direc- 
tion, y//j,2,o, and in the imaginary direction, JjMyfl, see 

Eq. (25]), both scale with V / N(R} 2 . This means that the 
typical decay rate is (A)' = N(R), whereas the typical 
oscillation frequency is of the order of y/~N (R). For large 
TV it will thus be difficult to observe the oscillations. 
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FIG. 5: (Color online) Fraction /r of nonzero eigenvalues that 
are real, as a function of TV for the EGRE. The solid cir- 
cles denote numerical values obtained for n r realizations with 
n r N = 4 x 10 7 for TV < 2000, n r N = 10 7 for TV = 5000, and 
n r N = 4 x 10 s for TV = 10000. The solid square represents 
the exact result / R = 1 for TV = 2. The dashed line denotes a 
power law foN~ a fitted to the two points for TV = 2000 and 
TV = 5000. The solid line is the exact result for the GinOE, 
Eq. (|29[) . The crosses denote numerical results for ensembles 
of rate matrices with Gaussian instead of exponential distri- 
bution of rates Aij, i 7^ j. 
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FIG. 7: (Color online) Scaled distribution function of complex 
eigenvalues A of shifted general rate matrices A of dimension 
(a) TV = 100 and (b) TV = 2000. Specifically, two-dimensional 
histograms with 500 x 500 bins were populated for n r matrices 
randomly chosen from the EGRE, where n r TV = 4 x 10 7 . Note 
the different scales of the axes. 



FIG. 6: (Color online) Scaled density of nonzero real eigenval- 
ues of shifted general rate matrices A. The curve for TV = 2 is 
exact. The curves for TV = 10, 100, 1000, 5000 are histograms 
with 500 bins for a total number n r N = 4 x 10 7 (10 7 ) of 
eigenvalues (including complex ones) for TV < 1000 (5000) for 
matrices randomly chosen from the EGRE. 



It is instructive to compare the distribution to the one 
for the GinOE. For the GinOE, the distribution function 
PC of complex eigenvalues for finite TV has been obtained 
by Edelman [3l| in terms of a finite sum of TV — 1 terms, 
which can be rewritten as a simple integral [33| . The dis- 
tribution function pc is found to contain a factor |ImA|, 
showing that the density goes to zero linearly for A ap- 
proaching the real axis. Complex eigenvalues are thus 
repelled by the real axis with a characteristic exponents 



of unity. Figures|3]and[7]clearly show that complex eigen- 
values are also repelled by the real axis for the EGRE. 
In Fig. [8] we plot the density of complex eigenvalues, pro- 
jected onto the real and imaginary axes, for TV = 100 
and TV = 2000. We observe that for the EGRE the com- 
plex eigenvalues are repelled by the real axis with the 
same exponent of unity. We note that the distribution 
of the real part of complex eigenvalues is distinct from 
both the distribution of real eigenvalues, Fig. O and the 
distribution of eigenvalues for the ESRE, Fig. [TJ 

For the GinOE, the scaled distribution approaches a 
uniform distribution on the unit disk in the c omp lex 
plane for TV — > 00. This was conjectured by Girko [28| for 
an arbitrary distribution of components with zero mean 
and proven by Bai [3(| • The EGRE result is clearly much 
more complicated. The histograms for various values of 
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FIG. 8: (Color online) Scaled density of the real part (solid 
lines) and the imaginary part (dashed lines) of complex eigen- 
values of shifted general rate matrices with N = 100 and 
N = 2000. The curves are projections of the data shown in 
Fig. [7] onto the real and imaginary axes. The inset shows the 
scaled density of the imaginary part around zero. 



N suggest that the distribution function pc does not be- 
come uniform in a bounded region for N — > oo, although 
it does appear to develop non-analyticities, which show 
up as high-contrast edges in Fig. EJb). 

We now return to the moments of the distribution func- 
tion p of all nonzero eigenvalues of A. The moments p mn , 
Eq. (j2"5|) , and the pseudomoments p m , Eq. (fM| . are re- 
lated. This is easily seen for p2- 

p 2 = (A 2 )' = ((Re A) 2 + 2i Re Aim A - (ImA) 2 ). (31) 

Since the second term vanishes, we obtain p2 = M2,o — 
/io,2- Now p2,o contains contributions from the real and 
the complex eigenvalues, while ^0,2 only depends on the 
complex eigenvalues. We can write 



A*2 = fa f4 + (1 - fa) M2.0 ^ (1 



e \ C 

JKJ Mo,2> 



(32) 



where the superscript K or C refers to the moments of 
the distributions of real and complex eigenvalues, respec- 
tively. In the limit of large N we know that /r — » and 
M2 — N(R) 2 . This means that the scaled distribution 
in the complex plane must be anisotropic: The width in 
the imaginary direction must be smaller by a value of the 
order of unity than in the real direction, unlike for the 
GinOE. This is seen in Fig. [7] 

For arbitrary even m, the relation reads 



i/2 



/'.;; ^ ^ ( 1) 1 ^ J /^m — n,n 

m 

/«A + (l-/t) E (-I)" 72 Q /£-„,„■ 

(33) 



n— even 



We recall that the p m for small to are known for all N, see 
Table HV1 For large N, we have the asymptotically exact 



expression (|C8|) . which can be written as p m = (m — 
1)!! iV m / 2 (R) m . Hence, we find asymptotically exact sum 
rules for all even orders to. 

To end this section, we again consider the slowest pro- 
cess. The dynamics at late times is typically governed 
by the eigenvalue Ai with the largest (smallest in magni- 
tude) real part. In Fig. [5] we show the mean of the real 
part Re Ai and of the magnitude of the imaginary part, 
|Im Ai I for random matrices from the EGRE, as functions 
of N. The behavior of the real part, i.e., the rate, is very 
similar to the ESRE. Again, the slowest rate is consistent 
with the mean and width of the eigenvalue distribution 
p(A). The typical imaginary part of Ai, i.e., the oscil- 
lation frequency, decreases for large N, mainly because 
the probability of Ai being real increases. While the frac- 
tion of real eigenvalues approaches zero for large N, the 
eigenvalue with the largest real part becomes more likely 
to be real. 
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FIG. 9: (Color online) Typical real and imaginary parts of the 
eigenvalue belonging to the slowest non-stationary process, as 
functions of N. The open circles denote the average smallest 
in magnitude real part of eigenvalues, (ReAi), of matrices 
from the EGRE. The data are numerical results for n r — 5000 
(1000) realizations for N < 1000 (N > 2000). The crosses 
denote the average magnitude of the imaginary part of the 
same eigenvalues, (|ImAi|), scaled xlOO. Error bars denoting 
the statistical errors are shown. The filled square denotes the 
exact result Ai = — 2{B) for N = 2. The dashed line denotes 
the mean of nonzero eigenvalues, —N{R). 



B. Eigenvalue correlations 

The eigenvalue density for the EGRE is quite differ- 
ent from the GinOE. Like for the ESRE, we again ask 
whether the eigenvalue correlations are also different. 
We consider the real and complex eigenvalues separately. 
The main effect of correlations between real and complex 
eigenvalues is seen in Fig. [5] The complex eigenvalues are 
repelled by the real axis with a characteristic exponent 
of unity. 
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Fi gure 1101 shows the distribution function pj\jN (AA) of 
separations of neighboring real eigenvalues. Note that the 
distribution is not rescaled with a power of N. The typ- 
ical separation of real eigenvalues depends only weakly 
on N for large N for the EGRE, whereas it scales with 
N- 1 ' 2 for the ESRE. This can be understood as fol- 
lows: The expected number of real eigenvalues of a ran- 
domly chosen matrix is Nfg. ~ iV 1_Q , while the width 
of their distribution scales with iV 1/,2 (i?). Consequently, 
the typical nearest-neighbor separation should scale with 
N a "^ 2 (R). Since a is close to 1/2, we obtain a weak 
dependence on N. The dependence on separation AA is 
again linear for small AA, though. Thus real eigenvalues 
repel each other with a characteristic exponent of unity, 
like for the GinOE 1331. 




4 6 
Al/<R> 

FIG. 10: (Color online) Distribution of nearest-neighbor sep- 
arations AA of nonzero real eigenvalues for the EGRE for 
N = 10, 100, 1000, 5000 for the same data sets as in Fig. [6] 
The axes are not rescaled with a power of N. 

In Figs. [TIT a) and (b), we plot the distribution function 
Pnn^'M °f complex differences of neighboring complex 
eigenvalues with positive imaginary part for TV = 20 and 
N = 2000. More specifically, for each eigenvalue A with 
positive imaginary part, we determine the eigenvalue A' 
with positive imaginary part that minimizes |A' — A|. We 
then collect the complex differences AA = A' — A of all 
such pairs in a two-dimensional histogram. The eigen- 
values with negative imaginary part just form a mirror 
image. Correlations between eigenvalues with positive 
and negative imaginary parts are dominated by their re- 
pulsion by the real axis and a (^-function from complex 
conjugate pairs and are not considered further. 

Since the fraction of complex eigenvalues approaches 
unity for TV — > oo, the number of complex eigenvalues 
of a chosen matrix scales with TV. The widths of the 
distribution in both the real and the imaginary direction 
scale with y/N, see Fig. The typical nearest-neighbor 
distance should thus approach a constant for large N. 
This is indeed seen in Fig. [TTJ 

We observe that the distribution of differences becomes 
rotationally symmetric for large N. This is perhaps sur- 
prising since the distribution of the eigenvalues them- 
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FIG. 11: (Color online) Distribution function of complex dif- 
ferences AA of neighboring eigenvalues with positive imagi- 
nary part for the EGRE for (a) N = 20 and (b) N = 2000. 



selves is far from symmetric, see Fig. [7J Also, small dif- 
ferences are suppressed, i.e., the eigenvalues repel each 
other. To find the characteristic exponent, we plot the 
distribution of the magnitudes |AA| = |A' — A| of differ- 
ences of neighboring eigenvalues in Fig. [T3J We observe 
that the distribution behaves like |AA| 3 for small |AA|. 
Together with the rotational symmetry this implies that 
the two-dimensional distribution in the complex plane, 
Fig. fTTT b). approaches zero like |AA| 2 . The exponent 
of two is the same as for the GinOE [33J . We conclude 
that the constraint ([5]) and the exponential distribution 
of rates in the EGRE do not change the repulsion of 
neighboring eigenvalues compared to the GinOE, while 
the eigenvalue density is very different. The origin of this 
is likely the same as to the ESRE: The correlations are 
governed by "local" properties of the joint distribution 
function of eigenvalues, which are not strongly affected 
by the constraint. 
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FIG. 12: (Color online) Distribution function of nearest- 
neighbor distances |AA| for various values of N. The dashed 
curve shows a power law oc |AA| 3 . 



IV. CONCLUSIONS 

We have applied RMT to the transition-rate matrix A, 
i.e., the matrix of coefficients in the Pauli master equa- 
tion ([3]) . This allows us to obtain statistical properties of 
the spectrum, in analogy to RMT for Hamiltonians. For 
the master equation, the eigenvalues describe the decay, 
and, in the case of complex eigenvalues, the superimposed 
oscillations, of probability eigenvectors. 

The resulting random-matrix ensembles are different 
from the standard ensembles for Hamiltonians, since A 
is real but in general not symmetric and since the conser- 
vation of probability imposes the constraint Aij = 
for all j, Eq. ([5]). Although this constraint represents only 
N conditions for of the order of N 2 matrix components, 
its consequences persist for large N. 

A further difference to the standard ensembles is that 
the off-diagonal components of the rate matrix represent 
rates and thus must be non-negative. We have assumed 
an exponential distribution. The results in the large-iV 
limit are found to be independent of the distribution of 
rates, though. 

We have considered both symmetric and general, 
asymmetric rate matrices. The first case corresponds to 
systems where the rates for transitions from any state i to 
any other state j and from j to i are identical. In the sec- 
ond case, these rates are assumed to be independent. In 
both cases, all nonzero eigenvalues form a narrow distri- 
bution of width proportional to y/~N around their mean, 
—N(R), where (R) is the average transition rate. Thus 
for not too small N, nearly all deviations from the sta- 
tionary state decay on the same time scale 1/N(R). For 
both cases, we have found that the slowest non-stationary 
state, which dominates the dynamics at late times, typ- 
ically also decays on the same time scale. We have de- 
rived exact expressions for the expectation values of m-th 
powers of the nonzero eigenvalues, for small m, for both 
cases. 



For symmetric rate matrices, the density of eigenval- 
ues has been studied numerically as a function of N and 
found to approach the same limiting form for N — > oo 
as obtained earlier for Gaussian and two- valued distribu- 
tions [22J, |23, l24j , but very different from the semi-circle 
law for the GOE [H QjJ. This difference is due to the 
constraint (JSJ. On the other hand, the correlations be- 
tween eigenvalues are dominated by a repulsion with a 
characteristic exponent of unity, as for the GOE. 

For general rate matrices, we have numerically studied 
the eigenvalue density in the complex plane as a func- 
tion of N. For large N, it approaches a non-trivial dis- 
tribution different from the disk found for the GinOE 
[H, |3(| • Interestingly, the fraction of nonzero eigenvalues 
that are real decays as N~ a with an anomalous exponent 
a w 0.460, unlike for the GinOE, where a = 1/2. Thus 
the fraction of eigenvectors describing purely exponen- 
tially decaying deviations from the stationary state scales 
with a nontrivial power of the number of possible states. 
Both the non-trivial distribution and the anomalous scal- 
ing for large N are due to the constraint (j5]) . The density 
of real eigenvalues is also different from the GinOE. We 
have obtained simple analytical results for the expecta- 
tion values ((A - (A)')™)' = (m - 1)!! (N(SR 2 )) m / 2 of all 
even powers of shifted nonzero eigenvalues in the limit 
of large N. Interestingly, they agree with the central 
moments of a real Gaussian distribution. The central 
moments of the eigenvalue density in the complex plane 
are shown to satisfy exact sum rules involving these ex- 
pectation values. 

Correlations between eigenvalues are found to agree 
with the GinOE: Real eigenvalues repel each other with 
an exponent of unity, complex eigenvalues are repelled 
by the real axis with an exponent of unity and by each 
other with an exponent of two. 

In view of the power of RMT for Hamiltonians, we 
hope that this approach will also benefit our understand- 
ing of complex stochastic processes. Comparisons with 
real processes are now called for. 



APPENDIX A: LARGE- iV LIMIT FOR 
SYMMETRIC RATE MATRICES 

In the limit of large N, the density of eigenvalues A of 
A only depends on the second moment (SR 2 ) of the dis- 
tribution of components Aij, i ^ j, for any distribution 
function of A^, as long as all its central moments exist. 
In this appendix, we sketch the proof of this statement. 

The eigenvalue density is given by Eq. (p~8|) . In the 
expansion of the geometric series for the resolvent [26[ , 

the n = term is independent of the distribution of Ay , 
while the n = 1 term vanishes. Since X^(^ n )y = f° r 
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n > 1 we can write 



= J - E E 

1 00 1 



n=2 fcijfes,-" 



(A2) 



We now introduce a diagrammatic representation for the 
expectation values (A m )ij, i 7^ j: 

• < • s ^ = 0, (A3) 



< >< < * = E< i2 )y = E E^ ifc i)> 

(A4) 

m < >< < X < — • = ^ (I 3 ) y - etc. (A5) 



Here, an arrow represents a factor of A, a vertex (filled 
circle or cross) represents a matrix index, and all indices 
are summed over 1, . . . , TV, subject to the constraint that 
indices corresponding to filled circles are distinct. Ver- 
tices drawn as crosses do not imply any constraint. 

In Eq. (|A2[> . we now decompose the sums over indices 
into terms with equal and distinct indices. For equal in- 
dices we attach the arrows to the same filled-circle vertex, 
whereas distinct indices arc denoted by distinct filled- 
circle vertices. For example, 

E (^ 2 >« = * < >< < * 



• < — • < — • + Q-< — • + • < — Q 



(A6) 



The constraint Ajj — — Ylj-fj Aij assumes the form 

0. = 



o 
a 



(A7) 



where the open circle denotes an index that is differ- 
ent from the one connected to it but not otherwise con- 
strained. Applying this rule to all terms, we obtain open- 
circle vertices, which we dispose of by again distinguish- 
ing between equal and distinct indices. For example, 



E 



. . . - £\ 

• < — • < — • — • < — • — • < — • 



• < 



We have achieved that factors of A with two equal indices 
are no longer present and that all indices to be summed 
over are distinct. 

Since different off-diagonal components Aij are inde- 
pendent, except for Aji — A^, the expectation value of 
each term decays into a product of expectation values of 
powers of components, (SR m ) = ((Aij)" 1 ). The corre- 
sponding diagrams are of the forms 



0, 



(A9) 
(A10) 



= = 

= (SR 3 ) etc. (All) 





Finally, any term containing m vertices obtains a factor 
N(N~l)(N-2) ■ ■ ■ (JV-m+1) from the sum over distinct 
indices. In the limit of large N this becomes N m . 

We conclude that at any order n > 2 in Eq. (|A2[) . 
the largest terms for large N are the non-vanishing ones 
with the maximum number of vertices. Note that the 
diagrams generated by this procedure are always con- 
nected. Diagrams containing single arrows connecting 
two vertices vanish because of Eq. (|A9[) . For even n, 
the maximum number of vertices is n/2 + 1, which is 
obtained if all connections are double arrows. In this 
case the contribution is proportional to N n ^ 2+1 (SR 2 ) n ^ 2 . 
The next smaller terms have two triple arrows and con- 
tribute cx N n / 2 (5R 2 ) n / 2 ~ 3 (5R 3 ) 2 . For odd n, the largest 
terms have one triple arrow and all other connections 
are double arrows. Their contribution is proportional to 
N^ 2+1 / 2 (SR 2 )^ 2 -^ 2 (SR 3 ). 

Since Eq. (fT8|) contains an explicit factor of 1/N, 
the leading contributions to the density scale as TV"/ 2 
(TV"/ 2 " 1 / 2 ) for even (odd) n. If we rescale the density so 
that the width approaches a constant, the odd terms in 
the expansion (|A2[) vanish like TV -1 / 2 , showing that the 
rescaled density approaches an even function. Further- 
more, the leading even terms only depend on the second 
moment (SR 2 ), which is what we set out to prove. 

Rewriting Eq. (|A1[) in terms of the moments fi n , 



(G(z)) 



y n+l 



(A12) 



(A8) 



we see that the terms of order n contribute exclusively 
to the moment fi n . The result proved here is consistent 
with the calculated moments in Table |TJ 



APPENDIX B: LARGE- iV LIMIT FOR GENERAL 
RATE MATRICES 



For ensembles of general, asymmetric rate matrices, it 
is also true that the density of eigenvalues only depends 
on the second moment (SR 2 ) for large N. We here sketch 
the proof of this assertion. 
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The distribution of eigenvalues in the complex plane is 
given by Eqs. p5 ]H [2g ]) . We define 



j(r};z,z*) = Tr 



2N 




Z, Z 



(Bl) 



so that p a \\(x,y) = (1/irN) dg(0; z, z*)/dz* and expand 
the resolvent, 



Ddd 



X (I T - 2*7) 



(B2) 



Expanding the products, we obtain a linear combination 
of expressions of the form Tr (• • • A T ■ ■ ■ A - ■ ■) contain- 
ing any number of factors A T and A in any order. Now 
the arguments of App. \X\ go through with few changes. 
We can group the terms according to the total order 
to of A and A T . The term of order zero is indepen- 
dent of the distribution of Aij . The terms of first order 
are Tr (A) — Tr (A T ) = 0. In all other terms we can 
use cyclic permutation under the trace and the identity 
Tr B T = Tr B to make sure that a factor A and not A T is 
appearing first under the trace. We can then use Eq. ([5]) 
to write Tr (A ■••) = - Ey&M ~ '>«■ 

Now we can apply the diagrammatics of App. [X] 
(A T )ij = Aji is drawn as an arrow pointing in the op- 
posite direction. In the evaluation of expectation values 
corresponding to Eqs. (j AlOjl . (|A1 1|) we have to take into 
account that A^ and Aji are now independent so that 
we instead have 



o 

o=°< 



0, 

(5R 2 ), 





(SR 3 ), 




etc. 



(B3) 
(B4) 

(B5) 
(B6) 

(B7) 



We note that all terms of the same order to in Eq. (|B2|) 
have the same sign and thus cannot cancel. We thus find 
that to any order m the leading terms in Eq. (|B2j) for 
large N have the same form as for symmetric matrices. 
In particular, for even to the leading term in the den- 
sity pix, y) scales with N m f 2 {8R 2 ) m / 2 and the odd terms 
scale with a lower power of N. Finally, it is conceivable 
that taking the derivative of g(Q; z, z*) with respect to z* 
in order to obtain the density could remove the leading- 
N term. This is not the case, since for any even order 
to > 2 there is at least a contribution from to = n — 1 in 
Eq. (|B2[) . which is linear in z*. 



APPENDIX C: PSEUDOMOMENTS FOR THE 
EGRE 

In this appendix, we use the diagrammatics of App. [A] 
to calculate the pseudomoments 



Mm = (A™)' = 



1 



N — 1 



Tr (i r ' 



(CI) 



to > 2, to leading order for large TV for the EGRE. Ap- 
pendix [B] shows that for large N only the even pseudo- 
moments are relevant. We write 



N 



E < lw % 

l — E E (A kl A klk2 ---A km _ 1J ).(G2) 



It was shown in App. [B]that for large N the distribution 
of Aij only enters through its second moment (SR 2 ). We 
decompose all terms into a sum of contributions with 
equal or distinct indices, see Eq. (|A6[) . For each term, 
some or none of the indices in {i, k\, k%, . . . , fe m _i,j} are 
equal. Contributions for which two equal indices are sep- 
arated by other, distinct indices in this string, correspond 
to diagrams of the type 



and are of lower order in N. All remaining diagrams are 
of the form of chains leading from j to i with any number 
of single- vertex loops (A kk ) decorating the vertices. We 
call these single-vertex loops "leaves." 
Next, we prove 




= 



(C4) 



for N — > oo, where the left-most vertex in the first term 
carries I > 1 leaves, while the second vertex in the second 
term carries I — 1 > leaves. The shaded circle is an 
arbitrary diagram part. The proof proceeds as follows: 
Applying the rule (|A7|) . we obtain 



+ 



(C5) 



with the upper (lower) signs for even (odd) I. In the lead- 
ing large- N term, all connections must be of the form of 
two arrows pointing in the same direction, as in Eq. (|B4[) . 
This is only possible if we pair up the open-circle vertices 
among themselves, not with any vertices in the right- 
hand part of the diagrams. This requires / to be even. 
Furthermore, for the first diagram there are (I — 1)!! ways 
to partition I leaves into pairs. For the second diagram 
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there are I — 1 ways to pair one of the leaves with the left- 
most vertex and (/ — 3)!! ways to partition the remaining 
I — 2 leaves into pairs. With these factors we obtain 



(1-1 




= 0. (C6) 



All diagrams of leading order in N are of the form of 
one of the two diagrams in Eq. (|C4[) . Thus all diagrams 



cancel, except if only one of the two forms exists. This 
is only the case for 




= -(to- l)!!JV m / 2+1 {5R 2 ) m/2 , (C7) 

since its partner would contain only a single vertex, which 
is excluded by % ^ j . 

With the prefactor from Eq. (IC2|) . we obtain 

fi m = (to - 1)!! TV" 1 / 2 {5R 2 ) m/2 . (C8) 
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